﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "power")]
    public static unsafe double power(int n, IntPtr a_ptr, double eps, IntPtr v_ptr, int inter)
    {
        double* a = (double*)a_ptr.ToPointer();
        double* v = (double*)v_ptr.ToPointer();
        return power(n, a, eps, v, inter);
    }

    /// <summary>
    /// 乘幂法 求实矩阵的特征向量
    /// </summary>
    /// <param name="n"></param>
    /// <param name="a">a[n][n]实矩阵</param>
    /// <param name="eps">精度要求</param>
    /// <param name="v">v[n]特征向量</param>
    /// <param name="inter">最多迭代次数，1000</param>
    /// <returns>函数返回绝对值最大的特征值。在本函数程序返回时将显示迭代次数。本程序最多迭代1000次。</returns>
    public static unsafe double power(int n, double* a, double eps, double* v, int inter = 1000)
    {
        int i, j, k, flag = 1, iteration;
        double lambda, sum, z = 0.0, err, t = 0.0, d, f = 0.0;
        double* u = stackalloc double[n];

        iteration = 0;
        do
        {
            iteration++;
            // 计算u=Av
            for (i = 0; i < n; i++)
            {
                sum = 0.0;
                for (j = 0; j < n; j++)
                {
                    sum = sum + a[i * n + j] * v[j];
                }
                u[i] = sum;
            }
            d = 0.0;
            // 计算向量的范数			
            for (k = 0; k < n; k++)
            {
                d = d + u[k] * u[k];
            }
            d = Math.Sqrt(d);
            for (i = 0; i < n; i++)
            {
                v[i] = u[i] / d;
            }
            if (iteration > 1)
            {
                err = Math.Abs((d - t) / d);
                f = 1;
                if (v[0] * z < 0) f = -1;
                if (err < eps) { flag = 0; }
            }
            if (flag == 1)
            {
                t = d; z = v[0];
            }
            if (iteration >= inter)
            {
                flag = 0;
            }
        } while (flag == 1);

        lambda = f * d;
        return (lambda);
    }

#if __DRIVE_CODE__
      int main()
        {
          int i;
          double a1[3][3] = {{0,1,1.5},{-5,-0.5,1},{-1,2,3.5}};
          double a2[3][3] = {{-5,1,5},{1,0,0},{0,1,0}};
          double v[3] = {0,0,1};
          double lambda;
          lambda = power(3, &a1[0][0], 0.0000001, v);
          cout <<"绝对值最大的特征值 lambda1 = " <<lambda <<endl;
          for (i=0; i<3; i++)
                cout <<"v(" <<i <<")= " <<v[i] <<endl;
          cout <<endl;
          lambda = power(3, &a2[0][0], 0.0000001, v);
          cout <<"绝对值最大的特征值 lambda2 = " <<lambda <<endl;
          for (i=0; i<3; i++)
                cout <<"v(" <<i <<")= " <<v[i] <<endl;
          return 0;
        }
#endif
#if __DRIVE_CODE__
      int main()
        {
          int i;
          double a[3][3] = {{22,-121,400},{1,0,0},{0,1,0}};
          double v[3] = {0,0,1};
          double lambda;
          lambda = power(3, &a[0][0], 0.0000001, v);
          cout <<"绝对值最大的特征值 lambda = " <<lambda <<endl;
          for (i=0; i<3; i++)
                cout <<"v(" <<i <<")= " <<v[i] <<endl;
          cout <<endl;
          return 0;
        }
#endif
}